### Replication code for ################################################
### "The Strength of Party Labels in Brazil: Evidence from Survey Experiments"

### Requires the following files:
### 	replicfunctions.R (with the functions used here)
### 	_experimentvalid.RData (BEPS data)
### 	_experimentBEPSdata.RData (BEPS data)

### Make changes to paths in lines 25 and 26!!!!!!!!!!!!!!!

### Figures in the paper as saved as pdf files, as follows:
###		FIGURE 3A: fig-barplot13-combined.pdf
###		FIGURE 3B: fig-barplot45-combined.pdf
###		FIGURE 4a: fig-barplot13inout.pdf
###		FIGURE 4B: fig-barplot45inout.pdf
###		FIGURE 5A: fig-barplot0-combined-with13.pdf
###		FIGURE 5B: fig-barplot0-combined-with13.pdf

library(xtable)
library(car)
library(erer)
rm(list=ls(all=TRUE))

### Here, user much replace path to folders in local system
setwd("~/Dropbox/Data/Paper-PTExperiment") #folder where datafiles are
figure.folder <- "~/Dropbox/PT-Papers/Experiment paper/_Figures_Experiment_Paper/" 

source("__replicfunctions.R")
load("_experimentvalid.RData") 
load("_experimentBEPSdata.RData") # load BEPS data
Qs <- c("QVenezuela","QPresal","QEmpresas","QInvestimento","QSalario")

### TABLE 1 in the paper ###############
### Table the N by treatment:  #########
print(xtable(table(d$PP,d$T)[,c(1,4,2,3)]))
print(xtable(table(dbeps$PP,dbeps$T)))

### Analyze balance
if(1==2){
tmp <- na.omit(subset(d,select=c(TT,CEPRegion,R,Sexo,PP,Diff,Idade,iA)))
tmp$R <- car::recode(tmp$R,"c('08-15 SM','15+ SM')='08+ SM'")
library(Matching)                     
bal <- MatchBalance(I(TT=="_ctrl") ~ as.character(R)+CEPRegion+Sexo + PP+ Diff + Idade + iA,
             nboots=1000,data=tmp)
tt <- t(sapply(bal[[1]],extract))
rownames(tt) <- c(levels(tmp$R)[-1],levels(factor(tmp$CEPRegion))[-1],"Gender",levels(tmp$PP)[-1],levels(tmp$Diff)[-1],"Age","Activism")
print(xtable(tt,digits=2))

library(RItools)
tmp$TT <- tmp$TT=="_somecue"
thebal<-xBalance(TT~CEPRegion+R+Sexo+PP+Diff+Idade+iA,data=tmp,report="all")
plot(thebal,thecols=1,legend=F)
cat("Bowler and Hanson omnibus test\n");print(thebal$overall)
}


#### ESTIMATE MODELS WITH SEPPARATE TREATMENT CONDITIONS ####################
#### exp.eachtreat.fb() runs the models treating each treatment condition sepparately
#### effects.cue() computes "effects"
#### plot.effects() and plot.cue() makes plots
allTQs <- lapply(Qs,function(x){exp.eachtreat(x)}) ## Run 4 cue model for all variables
	## Extract each model from allTQs
	allTQs.pt.ctrl <- lapply(allTQs,function(x){x[["reg13ctrl"]]}) 
	allTQs.pt.noctrl <- lapply(allTQs,function(x){x[["reg13"]]})
	allTQs.psdb.noctrl <- lapply(allTQs,function(x){x[["reg45"]]})
	allTQs.psdb.ctrl <- lapply(allTQs,function(x){x[["reg45ctrl"]]})
	allTQs.pt.noctrl.bootse <- lapply(allTQs,function(x){x[["reg13boot"]]})
	allTQs.psdb.noctrl.bootse <- lapply(allTQs,function(x){x[["reg45boot"]]})
	allTQs.pt.ctrl.bootse <- lapply(allTQs,function(x){x[["reg13ctrlboot"]]})
	allTQs.psdb.ctrl.bootse <- lapply(allTQs,function(x){x[["reg45ctrlboot"]]})
names(allTQs.pt.noctrl)<-names(allTQs.psdb.noctrl)<- 
	names(allTQs.pt.ctrl)<-names(allTQs.psdb.ctrl)<-
	names(allTQs.pt.noctrl.bootse)<-names(allTQs.psdb.noctrl.bootse)<-
	names(allTQs.pt.ctrl.bootse)<-names(allTQs.psdb.ctrl.bootse )<-Qs

###### POOLING TREATMENTS TOGETHER  ###################################
#### exp.treatnot estimates the basic models with pooled cues
#### effects.pooled converts raw regression into "effects"
#### plot.pooled and plot.effects produces figures
allQs <- lapply(Qs,function(x){exp.sometreat(x)})#takes a while because of bootstraping
	allQs.tables <- lapply(allQs,function(x){x[["tables"]]})#extracts tables
	allQs.pt.ctrl <- lapply(allQs,function(x){x[["reg13ctrl"]]})
	allQs.pt.noctrl <- lapply(allQs,function(x){x[["reg13"]]})
	allQs.psdb.ctrl <- lapply(allQs,function(x){x[["reg45ctrl"]]})
	allQs.psdb.noctrl <- lapply(allQs,function(x){x[["reg45"]]})
	allQs.pt.noctrl.bootse <- lapply(allQs,function(x){x[["reg13boot"]]})#
	allQs.psdb.noctrl.bootse <- lapply(allQs,function(x){x[["reg45boot"]]})#
	allQs.pt.ctrl.bootse <- lapply(allQs,function(x){x[["reg13ctrlboot"]]})#
	allQs.psdb.ctrl.bootse <- lapply(allQs,function(x){x[["reg45ctrlboot"]]})#
	names(allQs.pt.noctrl)<-names(allQs.psdb.noctrl)<- 
	names(allQs.pt.ctrl)<-names(allQs.psdb.ctrl)<-
	names(allQs.pt.noctrl.bootse)<-names(allQs.psdb.noctrl.bootse )<-
	names(allQs.pt.ctrl.bootse)<-names(allQs.psdb.ctrl.bootse )<-Qs


#### PLOT AND ANALYZE RESULTS ##################################

	#Plot Effects on agreement with PT 
	#Figure in original submission, not in final version of the paper 
	plot.effects(effects.cue(allTQs.pt.noctrl,
							boot.se=allTQs.pt.noctrl.bootse))

	#Barplot presentation of each treatment
	#(version in paper includes pooled cues, below):
	#Results are almost identical with or without controls
	#No controls (pt and psdb)
	plot.cue(effects.cue(allTQs.pt.noctrl),13,cue="dbl")
	plot.cue(effects.cue(allTQs.pt.noctrl,
				boot.se=allTQs.pt.noctrl.bootse),13,cue="dbl")
	plot.cue(effects.cue(allTQs.pt.noctrl),13,cue="pt")
	plot.cue(effects.cue(allTQs.pt.noctrl),13,cue="psdb")
	plot.cue(effects.cue(allTQs.psdb.noctrl),45,cue="dbl")
	plot.cue(effects.cue(allTQs.psdb.noctrl,
				boot.se=allTQs.psdb.noctrl.bootse),45,cue="dbl")
	plot.cue(effects.cue(allTQs.psdb.noctrl),45,cue="psdb")
	plot.cue(effects.cue(allTQs.psdb.noctrl),45,cue="pt")
	#With controls (pt and psdb)
	plot.cue(effects.cue(allTQs.pt.ctrl,
				boot.se=allTQs.pt.ctrl.bootse),13,cue="dbl")
	plot.cue(effects.cue(allTQs.pt.ctrl),13,cue="pt")
	plot.cue(effects.cue(allTQs.pt.ctrl),13,cue="psdb")
	plot.cue(effects.cue(allTQs.psdb.ctrl),45,cue="dbl")
	plot.cue(effects.cue(allTQs.psdb.ctrl),45,cue="psdb")
	plot.cue(effects.cue(allTQs.psdb.ctrl),45,cue="pt")
	#Other partisans and non-partisans (double cue only)
	plot.cue(effects.cue(allTQs.pt.ctrl),0,cue="dbl") #for non partisans
	plot.cue(effects.cue(allTQs.pt.ctrl,
				boot.se=allTQs.pt.ctrl.bootse),66,cue="dbl") #for "others"
	
	### Effects of cues on agreement with PT (not in final paper)
	plot.effects.pooled(effects.pooled(allQs.pt.noctrl))
	plot.effects.pooled(effects.pooled(allQs.pt.noctrl),ommit.others=T)

	## PLOT barplots (plots actually in paper follow below)
	plot.cue(effects.pooled(allQs.pt.noctrl),13,cue="some") 
	plot.cue(effects.pooled(allQs.pt.noctrl,
				boot.se=allQs.pt.noctrl.bootse),13,cue="some")
	plot.cue(effects.pooled(allQs.psdb.noctrl),45,cue="some") 
	plot.cue(effects.pooled(allQs.psdb.noctrl,
				boot.se=allQs.psdb.noctrl.bootse),45,cue="some")#only for no ctrl
	plot.cue(effects.pooled(allQs.pt.ctrl,
				boot.se=allQs.pt.ctrl.bootse),13,cue="some") 
	plot.cue(effects.pooled(allQs.psdb.ctrl,
				boot.se=allQs.psdb.ctrl.bootse),45,cue="some") 
	
	## Barplot combining DOUBLE with SOME cue 
	## Without BEPS
	pdf(paste(figure.folder,"fig-barplot13.pdf",sep=""),width=8,height=8)
	plot.cue(effects.cue(allTQs.pt.noctrl,
				boot.se=allTQs.pt.noctrl.bootse),13,cue="dbl",
				x2=effects.pooled(allQs.pt.noctrl,
				boot.se=allQs.pt.noctrl.bootse))
	dev.off()
	pdf(paste(figure.folder,"fig-barplot45.pdf",sep=""),width=8,height=8)
	plot.cue(effects.cue(allTQs.psdb.noctrl,
				boot.se=allTQs.psdb.noctrl.bootse),45,cue="dbl",
				x2=effects.pooled(allQs.psdb.noctrl,
				boot.se=allQs.psdb.noctrl.bootse))
	dev.off()
	pdf(paste(figure.folder,"fig-barplot0.pdf",sep=""),width=8,height=8)
	plot.cue(effects.cue(allTQs.pt.noctrl,
				boot.se=allTQs.pt.noctrl.bootse),0,cue="dbl",
				x2=effects.pooled(allQs.pt.noctrl,
				boot.se=allQs.pt.noctrl.bootse))
	dev.off()

	## Barplot combining In and OUT (FIGURE 4 in PAPER)
	pdf(paste(figure.folder,"fig-barplot13inout.pdf",sep=""),width=8,height=8)
	plot.cue.H(effects.cue(allTQs.pt.noctrl,
				boot.se=allTQs.pt.noctrl.bootse),13,cue="pt",
				x2=effects.cue(allTQs.pt.noctrl,
				boot.se=allTQs.pt.noctrl.bootse),cue2="psdb")
	dev.off()
	pdf(paste(figure.folder,"fig-barplot45inout.pdf",sep=""),width=8,height=8)
	plot.cue.H(effects.cue(allTQs.psdb.noctrl,
				boot.se=allTQs.psdb.noctrl.bootse),45,cue="psdb",
				x2=effects.cue(allTQs.psdb.noctrl,
				boot.se=allTQs.psdb.noctrl.bootse),cue2="pt")
	dev.off()

##### Head-on test of party cues in FIGURE 4 #######
	## Are PT cues different than PSDB cues? No.
	## This was requested by last reviewer, and is mentioned in text:
	#Difference of PT and PSDB cues effects on PTista agreement with the PT
	tmppt <- sapply(allTQs.pt.noctrl.bootse,test.diff.effects,pty=13)  	
	#Difference of PT and PSDB cues effects on PSDBista agreement with the PSDB
	tmppsdb <- sapply(allTQs.psdb.noctrl.bootse,test.diff.effects,pty=45)  
	colnames(tmppt) <- colnames(tmppsdb) <- Qs
	
cat("The absolute difference in the effects of PT and PSDB cues reported in both panels in Figure ref{figbarinout} is never larger than ",round(max(abs(c(tmppt["diff",],tmppsdb["diff",]))),2),", the mean absolute difference is much smaller (",round(mean(abs(c(tmppt["diff",],tmppsdb["diff",]))),2),"), and p-values for these differences computed using bootstrapped standard errors are are always greater than ",round(min(abs(c(tmppt["pval",],tmppsdb["pval",]))),3),".",sep="")
	
####### Add BEPS results to plots for paper #####
#### The technical problem here is that 'some' in BEPS is really 'dbl' in FB (for historical reasons)
#### So BEPS is anlayzed using the "pooled" functions (for a single treatment)
#### But has to be reported using the "cue" functions (for multiple treatments, as reviewer poitned out)
#### This is achieved by manually adding BEPS results to effects.cue objects
#### The result is cluncky (below) but avoids extensive re-writing of all funcitons

	#Estimate BEPS results and produce "effects" object to be merged with FB effects
	allQs.beps<-exp.sometreat("Qempresas",the.data=dbeps) #BEPS (labels is "some")
	beps.pt.noctrl <- effects.pooled(list(allQs.beps[['reg13']]),
						boot.se=list(allQs.beps[['reg13boot']]))
	beps.psdb.noctrl <- effects.pooled(list(allQs.beps[['reg45']]),
						boot.se=list(allQs.beps[['reg45boot']]))

	#Merge results for ploting for the PT
	comb.pt <- effects.cue(allTQs.pt.noctrl,boot.se=allTQs.pt.noctrl.bootse)
	comb.pt$baseline <- rbind(comb.pt$baseline,
						'QEmpresasBEPS'=c(beps.pt.noctrl$baseline)) #add baseline
	comb.pt$pes <- rbind(comb.pt$pes,
				'QEmpresasBEPS'=c(NA,NA,beps.pt.noctrl$pes[1],
				NA,NA,beps.pt.noctrl$pes[2],
				NA,NA,beps.pt.noctrl$pes[3],beps.pt.noctrl$pes[4]))
	comb.pt$ses <- rbind(comb.pt$ses,	
				'QEmpresasBEPS'=c(NA,NA,beps.pt.noctrl$ses[1],
				NA,NA,beps.pt.noctrl$ses[2],
				NA,NA,beps.pt.noctrl$ses[3],beps.pt.noctrl$ses[4]))
	comb.pt.pooled <- effects.pooled(allQs.pt.noctrl,
					boot.se=allQs.pt.noctrl.bootse)
	comb.pt.pooled$baseline <- rbind(comb.pt.pooled$baseline,'QEmpresasBEPS'=c(NA)) 
	comb.pt.pooled$pes <- rbind(comb.pt.pooled$pes,'QEmpresasBEPS'=c(NA)) 
	comb.pt.pooled$ses <- rbind(comb.pt.pooled$ses,'QEmpresasBEPS'=c(NA)) 

	#Merge results for plotting for the PSDB
	comb.psdb <- effects.cue(allTQs.psdb.noctrl,boot.se=allTQs.psdb.noctrl.bootse)
	comb.psdb$baseline <- rbind(comb.psdb$baseline,
				'QEmpresasBEPS'=c(beps.psdb.noctrl$baseline)) 
	comb.psdb$pes <- rbind(comb.psdb$pes,
				'QEmpresasBEPS'=c(NA,NA,beps.psdb.noctrl$pes[1],
				NA,NA,beps.psdb.noctrl$pes[2],
				NA,NA,beps.psdb.noctrl$pes[3],beps.psdb.noctrl$pes[4]))
	comb.psdb$ses <- rbind(comb.psdb$ses,	
				'QEmpresasBEPS'=c(NA,NA,beps.psdb.noctrl$ses[1],
				NA,NA,beps.psdb.noctrl$ses[2],
				NA,NA,beps.psdb.noctrl$ses[3],beps.psdb.noctrl$ses[4]))

	comb.psdb.pooled <- effects.pooled(allQs.psdb.noctrl,
					boot.se=allQs.psdb.noctrl.bootse)
	comb.psdb.pooled$baseline <- rbind(comb.psdb.pooled$baseline,'QEmpresasBEPS'=c(NA)) 
	comb.psdb.pooled$pes <- rbind(comb.psdb.pooled$pes,'QEmpresasBEPS'=c(NA)) 
	comb.psdb.pooled$ses <- rbind(comb.psdb.pooled$ses,'QEmpresasBEPS'=c(NA)) 

	### FIGURE 3 in paper #################
	#Control, Some, and Dble cue (function.H)
	pdf(paste(figure.folder,"fig-barplot13-combined.pdf",sep=""),width=8,height=8)
	plot.cue.H(comb.pt,x2=comb.pt.pooled,13,cue="dbl")
	dev.off()
	pdf(paste(figure.folder,"fig-barplot13-combined-dblonly.pdf",sep="")
			,width=8,height=8)
	plot.cue.H(comb.pt,13,cue="dbl")
	dev.off()
	
	pdf(paste(figure.folder,"fig-barplot45-combined.pdf",sep=""),width=8,height=8)
	plot.cue.H(comb.psdb,x2=comb.psdb.pooled,45,cue="dbl")
	dev.off()
	pdf(paste(figure.folder,"fig-barplot45-combined-dblonly.pdf",sep="")
		,width=8,height=8)
	plot.cue.H(comb.psdb,45,cue="dbl")
	dev.off()
	
	### FIGURE 5 in paper #################
	pdf(paste(figure.folder,"fig-barplot0-combined.pdf",sep=""),width=8,height=8)
	plot.cue.H(comb.pt,x2=comb.pt.pooled,0,cue="dbl")
	dev.off()
	pdf(paste(figure.folder,"fig-barplot0-combined-with13.pdf",sep=""),width=8,height=8)
				plot.cue.H(comb.pt,0,cue="pt",x2=comb.pt,cue2="dbl")
	dev.off()
	pdf(paste(figure.folder,"fig-barplot0-combined-with45.pdf",sep=""),width=8,height=8)
		plot.cue.H(comb.psdb,0,cue="psdb",x2=comb.psdb,cue2="dbl")
	dev.off()


#report average effects across questions with less than .6 baseline agreement
	cat("Considering only these low-agreement questions, the average double-cue treatment effect across all questions is actually smaller for petistas (",round(mean(comb.pt[['pes']][which(comb.pt[['baseline']][,"PP13"]<0.6),"PP13:T4dblcue"]),2),") than for PSDB supporters (",round(mean(comb.psdb[['pes']][which(comb.psdb[['baseline']][,"PP45"]<0.6),"PP45:T4dblcue"]),2),"). If we exclude the BEPS results, the average effects become almost identical --- ",round(mean(comb.pt[['pes']][which(comb.pt[['baseline']][-6,"PP13"]<0.6),"PP13:T4dblcue"]),2)," for supporters of the PT, and ",round(mean(comb.psdb[['pes']][which(comb.psdb[['baseline']][-6,"PP45"]<0.6),"PP45:T4dblcue"]),2)," for those of the PSDB. Our conclusion is that for similar levels of baseline agreement, the effects of the PSDB label on sympathizers of the PSDB are similar to that of the PT label on its own sympathizers.",sep="")


#### OTHER Parties (NOT INCLUDED IN FINAL PAPER) ######################################
### Estimate regressions unbundling "others"
### Takes a long time to run!
allTQs66 <- lapply(Qs,function(x){exp.eachtreat(x,partyvar="P")}) 
allQs66 <- lapply(Qs,function(x){exp.sometreat(x,partyvar="P")})

	#Extract results
	allQs66.pt.noctrl <- lapply(allQs66,function(x){x[["reg13"]]})
	allQs66.pt.noctrl.bootse <- lapply(allQs66,function(x){x[["reg13boot"]]})
	allQs66.pt.ctrl <- lapply(allQs66,function(x){x[["reg13ctrl"]]})
	allQs66.pt.ctrl.bootse <- lapply(allQs66,function(x){x[["reg13ctrlboot"]]})
	
	allTQs66.pt.noctrl <- lapply(allTQs66,function(x){x[["reg13"]]})
	allTQs66.pt.noctrl.bootse <- lapply(allTQs66,function(x){x[["reg13boot"]]})
	
	pdf(paste(figure.folder,"fig-barplot43-pooledonly.pdf",sep=""),width=8,height=8)
	plot.cue(effects.pooled.other(allQs66.pt.noctrl,
				boot.se=allQs66.pt.noctrl.bootse),43,cue="some")
	dev.off()
	pdf(paste(figure.folder,"fig-barplot15-pooledonly.pdf",sep=""),width=8,height=8)
	plot.cue(effects.pooled.other(allQs66.pt.noctrl,
				boot.se=allQs66.pt.noctrl.bootse),15,cue="some")	
	dev.off()
	pdf(paste(figure.folder,"fig-barplot66-pooledonly.pdf",sep=""),width=8,height=8)
	plot.cue(effects.pooled.other(allQs66.pt.noctrl,
				boot.se=allQs66.pt.noctrl.bootse),66,cue="some")		
	dev.off()		
	
	pdf(paste(figure.folder,"fig-effects-other.pdf",sep=""),width=8,height=8)
	plot.effects.pooled(effects.pooled.other(allQs66.pt.ctrl,
							boot.se=allQs66.pt.ctrl.bootse),leg2=F)
	legend(x="bottom",
		legend=c("PMDB","Other","Non-Partisan","PV"),
		cex=1.2,pch=c(21,21,23,21),pt.bg=c(1,gray(0.5), gray(0.5),"white"),
		col=c(1,gray(0.5),gray(0.5),1),pt.cex=2,xpd=NA,
		inset=c(0,-0.23),bty="n",horiz=T)
	dev.off()
			
#### INTERNAL VALIDITY ISSUES  (in Supplemental Information)
#### Party ID given ordering of ID question   ########################
d$Tid <- factor(ifelse(d$Qconditionpid==1,"_1IDbeforeEXP","_2IDafterEXP"))
d$PPP <- d$PP!=0 #identifies with PT or PSDB

reg1 <- summary(lm(PPP~Tid,data=d)) #id with PT and PSDB given ID before or after
reg2 <- summary(lm(PPP~Tid,data=subset(d,TT=="_somecue"))) #id with PT and PSDB given ID before or after, having received cue
reg4 <- summary(lm(I(PP==13)~Tid,data=subset(d,TT=="_somecue"))) 
reg5 <- summary(lm(I(PP==45)~Tid,data=subset(d,TT=="_somecue"))) 

tmp <- data.frame(rbind(
reg1$coef[2,],
reg2$coef[2,],
reg4$coef[2,],
reg5$coef[2,]
))
print(xtable(tmp[,1:2],digits=3))
summary(lm(I(d$Diff!='0Nenhuma')~Tid,data=d))

### Cueing effects given position of ID question   ########################
old.d <- d  #keep copy of d, because function calls general (not local) d
d <- subset(old.d,Tid=="_1IDbeforeEXP")
noctrl.bef <- lapply(Qs,function(x){exp.treatnot(x)[[3]]})
effects.before <- plot.pooled.treatments(noctrl.bef)

d <- subset(old.d,Tid=="_2IDafterEXP")
noctrl.after <- lapply(Qs,function(x){exp.treatnot(x)[[3]]})
effects.after <-  plot.pooled.treatments(noctrl.after)

pes.before <- do.call("rbind",effects.before[['pes']])
pes.after <- do.call("rbind",effects.after[['pes']])

diff.befafter <- pes.after-pes.before
rownames(diff.befafter) <- Qs
diff.befafter <- rbind(diff.befafter, apply(diff.befafter,2,mean))
print(xtable(diff.befafter,digits=3))